!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine mpact (et, ro, mx, my, bx, by, en, az, ro_old, mx_old, my_old, dx, dy, dt, nx, ny)
!==================================================================================================================================|

   implicit none

   integer(4), intent(in) :: nx, ny
   real(8), intent(in) :: dx, dy, dt
   real(8), dimension(nx, ny), intent(in) :: ro, mx, my, ro_old, mx_old, my_old
   real(8), dimension(nx, ny), intent(inout) :: et, bx, by, en, az

   integer(4) :: i, j, ip1, jp1, im1, jm1
   real(8) :: dt_dx, dt_dy, dx2, dy2, dxdx, dydy, wx, wy, vxh, vyh, theta, eps
   real(8), dimension(nx, ny) :: bx_old, by_old
   real(8), dimension(nx, ny) :: vx, vy, az_pre, lim_x, lim_y, dw_x, dw_y
   real(8), dimension(nx, ny) :: fcor1_m_x, fcor1_p_x, fcor1_m_y, fcor1_p_y, fcor2_m_x, fcor2_p_x, fcor2_m_y, fcor2_p_y

!----------------------------------------------------------------------------------------------------------------------------------|
   dx2 = 2.0d0 * dx
   dy2 = 2.0d0 * dy
   dxdx = dx * dx
   dydy = dy * dy
   dt_dx = dt / dx
   dt_dy = dt / dy
   eps = 1.0d-12

   do j = 1, ny
   do i = 1, nx
      bx_old(i, j) = bx(i, j)
      by_old(i, j) = by(i, j)
      vx(i, j) = (mx(i, j) / ro(i, j) + mx_old(i, j) / ro_old(i, j)) * 0.5d0
      vy(i, j) = (my(i, j) / ro(i, j) + my_old(i, j) / ro_old(i, j)) * 0.5d0
   enddo
   enddo

!   do j = 2, ny - 1
!   do i = 2, nx - 1
!      az_pre(i, j) = az(i, j) - dt_dx2 * vx(i, j) * (az(i + 1, j) - az(i - 1, j)) -                                                &
!         dt_dy2 * vy(i, j) * (az(i, j + 1) - az(i, j - 1))
!   enddo
!   enddo

   do j = 2, ny - 1
   do i = 2, nx - 1
      dw_x(i, j) = az(i + 1, j) + az(i - 1, j) - 2.0d0 * az(i, j)
      dw_y(i, j) = az(i, j + 1) + az(i, j - 1) - 2.0d0 * az(i, j)
   enddo
   enddo

   do j = 2, ny - 1
   do i = 2, nx - 1
      if (vx(i, j) .ge. 0) then
         theta = dw_x(i - 1, j) / (dw_x(i, j) + eps)
      else
         theta = dw_x(i + 1, j) / (dw_x(i, j) + eps)
      endif
      lim_x(i, j) = max (0.0d0, min (min ((1.0d0 + theta) / 2.0d0, 2.0d0), 2.0d0 * theta))
!      lim_x(i, j) = max (0.0d0, min (1.0d0, theta))
!      if (abs (theta) .lt. 1.0d0) then
!         lim_x(i, j) = theta
!      else
!         lim_x(i, j) = 0.0d0
!      endif

      if (vy(i, j) .ge. 0) then
         theta = dw_y(i, j - 1) / (dw_y(i, j) + eps)
      else
         theta = dw_y(i, j + 1) / (dw_y(i, j) + eps)
      endif
      lim_y(i, j) = max (0.0d0, min (min ((1.0d0 + theta) / 2.0d0, 2.0d0), 2.0d0 * theta))
!      lim_y(i, j) = max (0.0d0, min (1.0d0, theta))
!      if (abs (theta) .lt. 1.0d0) then
!         lim_y(i, j) = theta
!      else
!         lim_y(i, j) = 0.0d0
!      endif
   enddo
   enddo

   do j = 1, ny - 1
      jp1 = j + 1
      do i = 1, nx - 1
         ip1 = i + 1

         wx = az(ip1, j) - az(i, j)
         wy = az(i, jp1) - az(i, j)
         vxh = 0.5d0 * (vx(i, j) + vx(ip1, j))
         vyh = 0.5d0 * (vy(i, j) + vy(i, jp1))

!         fcor1_m_x(i, j) = min (vx(i, j), 0.0d0) * wx
!         fcor1_p_x(i, j) = max (vx(ip1, j), 0.0d0) * wx
!         fcor1_m_y(i, j) = min (vy(i, j), 0.0d0) * wy
!         fcor1_p_y(i, j) = max (vy(i, jp1), 0.0d0) * wy

!         fcor2_m_x(i, j) = 0.5d0 * (abs (vx(i, j)) - dt_dx * vxh * vx(i, j)) * wx * lim_x(i, j)
!         fcor2_p_x(i, j) = 0.5d0 * (abs (vx(ip1, j)) - dt_dx * vxh * vx(ip1, j)) * wx * lim_x(ip1, j)
!         fcor2_m_y(i, j) = 0.5d0 * (abs (vy(i, j)) - dt_dy * vyh * vy(i, j)) * wy * lim_y(i, j)
!         fcor2_p_y(i, j) = 0.5d0 * (abs (vy(i, jp1)) - dt_dy * vyh * vy(i, jp1)) * wy * lim_y(i, jp1)

         fcor1_m_x(i, j) = min (vxh, 0.0d0) * wx
         fcor1_p_x(i, j) = max (vxh, 0.0d0) * wx
         fcor1_m_y(i, j) = min (vyh, 0.0d0) * wy
         fcor1_p_y(i, j) = max (vyh, 0.0d0) * wy

         fcor2_m_x(i, j) = 0.5d0 * abs (vxh) * (1.0d0 - dt_dx * abs (vxh)) * wx * lim_x(i, j)
         fcor2_p_x(i, j) = 0.5d0 * abs (vxh) * (1.0d0 - dt_dx * abs (vxh)) * wx * lim_x(ip1, j)
         fcor2_m_y(i, j) = 0.5d0 * abs (vyh) * (1.0d0 - dt_dy * abs (vyh)) * wy * lim_y(i, j)
         fcor2_p_y(i, j) = 0.5d0 * abs (vyh) * (1.0d0 - dt_dy * abs (vyh)) * wy * lim_y(i, jp1)
      enddo
   enddo

   do j = 2, ny
      jm1 = j - 1
      do i = 2, nx
         im1 = i - 1
         az_pre(i, j) = az(i, j) - dt_dx * (fcor1_m_x(i, j) + fcor1_p_x(im1, j)) - dt_dx * (fcor2_m_x(i, j) - fcor2_p_x(im1, j)) - &
                                   dt_dy * (fcor1_m_y(i, j) + fcor1_p_y(i, jm1)) - dt_dy * (fcor2_m_y(i, j) - fcor2_p_y(i, jm1)) + &
                                   dt * et(i, j) * (dw_x(i, j) / dxdx + dw_y(i, j) / dydy)
      enddo
   enddo

   az = az_pre

   do j = 2, ny - 1
   do i = 2, nx - 1
      bx(i, j) = (az(i, j + 1) - az(i, j - 1)) / dy2
      by(i, j) = (az(i - 1, j) - az(i + 1, j)) / dx2
   enddo
   enddo

   do j = 1, ny
   do i = 1, nx
      en(i, j) = en(i, j) + 0.5d0 * (bx(i, j) * bx(i, j) + by(i, j) * by(i, j)) -                                                  &
         0.5d0 * (bx_old(i, j) * bx_old(i, j) + by_old(i, j) * by_old(i, j))
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine mpact
